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Abstract 

We  devise  a  class  of  fast  wavelet  based  algorithms  for  linear  evolution  equations  whose 
coefficients  are  time  independent.  The  method  draws  on  the  work  of  Beylkin,  Coifman,  and 
Rokhlin  [1]  which  they  applied  to  general  Calderon-Zygmund  type  integral  operators.  We 
apply  a  modification  of  their  idea  to  linear  hyperbolic  and  parabolic  equations,  with  spatially 
varying  coefficients.  A  significant  speedup  over  standard  methods  is  obtained  when  applied 
to  hyperbolic  equations  in  one  space  dimension  and  parabolic  equations  in  multidimensions. 
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1.  Introduction 


During  the  last  few  years  a  number  of  fast  computational  algorithms  have  been  developed 
for  elliptic  problems.  These  are  techniques  for  which  the  number  of  arithmetic  operations 
needed  are  close  to  linear  as  a  function  of  the  number  of  unknowns.  Examples  of  algorithms 
of  such  complexity  are  multigrid  methods  and  the  so-called  fast  Poisson  solvers.  The  fast 
multipole  method  and  wavelet  based  methods  for  elliptic  problems  formulated  as  integral 
equations  also  belong  to  this  category  [8],  [1]. 

There  has  not  been  the  same  progress  for  hyperbolic  and  parabolic  methods.  In  general 
classical  numerical  techniques  for  these  problems  are  already  optimal. 

Consider  a  system  of  evolution  equations. 

dtu  + L{x,dx)u  =  f(x),  a:  6  fl  C  <  >  0, 

(1.1) 

u(x  0)  =  Uq{x), 

with  boundary  conditions,  where  Z,  is  a  differential  operator. 

An  explicit  discretization  of  this  problem  typically  takes  the  form, 

Xj  =  {jiAxi,...,jd^Xd) 

(1.2)  = 

=  Uq, 

u,F  ^  At  =  const.  |Ax|’'. 

The  vector  u"  contains  all  the  unknowns  u"  at  time  level  For  simplicity  we  shall  assume 
=  1 , 2, . . . ,  iV  in  all  dimensions  i/  =  1, . . . ,  d. 

The  matrix  A  is  {N'^  x  N‘^)  with  the  number  of  elements  ^  0  in  each  row  and  each 
column  bounded  by  a  constant.  Every  time  step  requires  0{N'^)  arithmetic  operations  and 
the  overall  complexity  for  a  time  interval  of  0(\)  is  of  the  same  order  as  the  number  of 
unknowns, 

There  are,  however,  some  fast  methods  based  on  the  analytic  form  of  the  solution  opera¬ 
tor.  In  [3]  the  multidimensional  heat  operator,  with  uq  and  /  both  zero,  but  with  inhomoge¬ 
neous  boundary  data  given  at  M  points,  was  treated.  There  the  closed  form  of  the  solution 
evaluated  at  M  points  at  time  level  N  v/as  obtained  in  0{NM)  rather  than  0{N^M‘^)  op¬ 
erations.  Also,  in  [4],  the  same  authors  obtained  an  algorithm  for  evaluating  the  sum  of 
N  Gaussians  at  M  arbitrarily  distributed  points  in  0{N  -f-  M)  operations.  So  far,  their 
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interesting  method  appears  to  need  an  explicit  analytic  representation  of  the  heat  kernel, 
effectively  ruling  out  variable  coefficient  problems. 

The  formula  (1.2)  has  a  simple  closed  form  solution 


(1.3) 


n— 1 

u"  =  A^Uq  +  Yi 

i/=0 


This  form  can  be  used  to  compute  the  solution  for  F  =  0,  in  log  n  steps,  (n  =  2*”,  m  in¬ 
teger;  here  and  throughout,  log  n  =  logj  n)  by  repeated  squaring  olA:  A,  A"^,  A"^,  /I®, . . . ,  A^”'. 

Unfortunately  the  later  squarings  involve  almost  dense  matrices  and  the  overall  complex¬ 
ity  is  0{N^‘^  log  N)  which  is  larger  than  that  using  (1.2)  directly. 

For  an  appropriate  representation  of  A  in  a  wavelet  basis  all  of  the  powers  A''  may  be 
approximated  by  spairse  matrices  and  the  algorithm  using  repeated  squaring  should  then  be 
advantageous. 

We  shall  consider  the  following  algorithms  for  the  computation  of  the  closed  form  solution 
(1.3)  of  the  inhomogeneous  problem  in  m  =  logn  steps. 


(1.4) 


B  :=  SAS-^ 

C:=l 

C  :=TRUNC{C  +  BC,€) 
B  :=  TRUNcIbB,£) 


(iterate  m  steps) 


u"  :=  S-^iBSu°  +  CSF). 

The  matrix  S  corresponds  to  a  fast  transform  of  wavelet  type  and  the  truncation  operator 
sets  elements  in  a  matrix  to  zero  if  their  absolute  value  is  below  a  given  threshold. 

(1.5)  A  =  TRUNC{A,e) :  (  ^  ^  ^  ^ 

It  is  easy  to  see  that  algorithm  (1.4)  is  equivalent  to  (1.3)  for  e  =  0.  This  is  not  so  for  £  >  0 
and  also  for  F  ^  0.  We  shall  however  show  that  it  is  possible  to  choose  e  small  enough  for 
the  result  of  (1.4)  to  be  arbitrarily  close  to  (1.3)  but  still  with  very  few  arithmetic  operations. 

For  a  fixed  predetermined  accuracy  level  the  computational  complexity  to  calculate  a  one 
dimensional  hyperbolic  equation  can  be  reduced  from  the  standard  O(iV^)  to  0{N{\og  N)^). 
The  extra  cost  per  time  step  is  minimal.  This  also  makes  it  possible,  as  a  curiosity,  to  use 
algorithms  which  are  unstable  in  the  traditional  sense. 

Our  technique  is  even  more  favorable  for  parabolic  problems.  A  d-dimensional  explicit 
calculation  with  standard  complexity  0(iV‘^'*'*)  may  be  reduced  to  0{N^{\ogN)^). 
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The  algorithm  (1.4)  can  be  extended  to  some  problems  with  time  dependent  data.  In  this 
case,  we  clearly  need  to  compress  the  information  in  the  data  such  that  not  all  the  0(iV'^+’') 
values  in,  e.g.  the  inhomogeneous  term  f{xj,tn)  are  needed. 

One  simple  but  important  application  of  this  type  is  from  optics  or  electro-magnetic 
scattering  with  a  time  periodic  source.  If  k  points  are  needed  to  resolve  one  time  period,  we 
can  group  k  time  steps  together 

(1.6a)  1*"+'=  =  -I-  2  A^Fn+k+j-i . 

i=o 

where 


(1.6b)  =  At/(<„). 

This  equation  is  now  of  the  type  (1.2)  with  time  step  kAt  and  with  inhomogeneous  term 

(1.6c) 


3=0 


In  sections  2  and  3  we  shall  discuss  the  analytical  properties  of  the  algorithm.  Numerical 
examples  are  presented  in  section  4. 


2.  Hyperbolic  Problems 

Consider  first  the  simple  one  dimensional  scalar  advection  equation, 

dtu  +adxU  =  0,  o  >  0 

u(x,0)  =  uo(x),  0  <  X  <  1. 


(2.1) 


The  functions  uq  eind  thus  u  are  assumed  to  be  1-periodic  in  x.  The  solution  of  (2.1)  is  given 
by: 


(2.2) 


u(x,t)  =  «o(x  —  at). 


The  different  rows  of  A*'  in  a  numerical  solution  of  (2.10)  will  represent  approximations  of 
the  Green’s  function  G  below. 


(2.3) 


/oo 

G{x,y,t)uo{y)dy, 

•OO 

/OO 

S{x-y-  at)uo{y)dy. 

•OO 
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Let  ifj  be  a  truncated  wavelet  expansion  of  a  5-function  with  an  orthonormal  set  of  compactly 
supported  wavelets, 

8{x)  ~  ipj{x)  =^Y^OLjk^jk{x) 

^jk{x)  =  2"^(2“-'x  -  A:  -t- 1) 

The  choices  of  ^l>{x)  will  be  discussed  below.  Assume  that  the  rows  of  A"  are  discrete  8- 
functions,  i.e.  just  one  element  is  nonzero  and  large.  For  each  level  j  =  1, 2, . . . ,  J  there  are 
only  a  finite  number  of  ajk  ^  0.  With  J  =  m  =  log  iV  there  is  only  log  N  of  all  Ojk  ^  0. 
Thus  each  row  in  J3,  (1.4),  has  logiV  elements,  bjk  ^  0.  The  matrix  is  also  a  transform 
of  an  idealized  matrix  A*'  and  will  have  N  log  N  elements  different  from  zero.  This  means 
that  each  iteration  step  in  the  algorithm  (1.4)  produces  0(A(logiV)^)  flops  when  F  =  0. 
We  have  assumed  that  calculations  are  only  carried  out  for  those  elements  which  are 
different  from  zero.  In  practice  a  slightly  larger  number  of  elements  needs  to  be  computed 
and  then  truncated.  This  corresponds  to  the  case  when  the  location  of  the  5-functions  is 
only  approximately  known.  Compare  the  wavelet  technique  for  Burgers’  equation  by  Maday, 
Perrier,  and  Ravel  [6]. 

Each  row  of  C,  (1.4),  is  a  transform  of  a  step  function, 

const.  0  <  X  <  at, 

0,  else 

This  function  can  also  be  represented  by  logiV  wavelets  and  thus  the  overall  cost  is 
0{N{log  N)^). 

In  numerical  computations  the  rows  of  A"  are  only  approximations  of  5-functions.  If  an 
upwind  scheme, 

=  u7-A(w7-u7_a), 
u°j  =  uo(a;,),  j  =  l,2,...N, 

X  =  aAt/Ax  <  1, 

is  used  A  will  have  the  form, 

1-A  0 

A  1-A  0 

0  A  1-A  0  ••• 

0  •  ■ '  0  A 

The  matrix  A"  will  have  Toeplitz  structure.  Each  row  is  still  an  approximation  of  a  5- 
function.  The  first  order  smoothing  effect  of  (2.4)  is  given  by  the  modified  equation,  [5], 

(2.5)  dtu  -f  adxU  =  {aAx/2)dlu. 
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Equation  (2.5)  is  parabolic  with  a  fundamental  solution  of  the  form, 


(2.6) 


G{x  —  y,t)  =  (27raAit)~2  exp(— (a:  —  y  —  ai)^/(2aAxf)). 


Compare  the  solution  formula  for  parabolic  problems  (3.2). 

Each  row  of  A''  is  thus  a  close  approximation  to  the  function  G{x  —  y,t)  above.  The 
computational  complexity  of  the  algorithm  (1.4)  depends  on  how  many  wavelets  are  needed 
to  represent  G{x  —  y,t)  as  &  function  of  x,(0  <  t  <  T)  with  a  given  accuracy. 

Higher  order  accurate  (say  order  2p-l)  dissipative  finite  difference  approximations  to  (2.1) 
are  usually  modelled  by  the  equation 


(2.7) 


/  d 

ut  +  aux  =  (-l)P'^^A:p(Ax)^’’ 


with  kp  >  6  >  0,6,  independent  of  Ax. 

The  fundamental  solution  for  this  parabolic  equation  is: 


0  =  ^  exp(i^(x  -  at)  -  fcp( Ax)^^  ^^^'’0)- 

The  key  estimate  we  shall  obtain  here  (and  which  we  certainly  do  not  claim  is  new)  is: 


(2.8) 


.m+l 


'  d 

-  G,(x  +  at,t) 


<C, 


uniformly  in  0  <  t  and  Ax  and  for  all  nonnegative  integers  m. 
Proof  of  2.8.  We  wish  to  bound 

27r  J—oo 


/  ^  \  n»+l 

J 

=  JL  r  fexp  ( _ [f  Ar" 

27ry_oo[  ^  V(t(Ax)2p-»A:p)V2pyJ  \d()  J 


d(. 


The  result  is  now  clear.  Also,  an  inspection  of  the  right  hand  side  of  the  above  shows  that 
C’m.p  can  be  chosen  to  be  arbitrarily  small  if  t(Ax)^'’  ^  is  large  enough. 


Remark  Rl.  Let  the  general  space  dependent  coefficient,  one  dimensional  system  of  hy¬ 
perbolic  equations 

Ut  +  A{x)ux  =  C{x)u, 
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where  w  is  an  ^  vector,  A  is  a  uniformly  diagonalizable  smooth  i  x  £  matrix,  with  all  real 
eigenvalues  A, (a;),  and  C{x)  is  smooth,  be  approximated  by  a  dissipative  finite  difference 
scheme  of  order  2p  —  1.  Typically,  its  model  equation  is  a  systems  version  of  (2.1) 

4)“ 

where  (— 1)p+^P(x,  is  a  2p  order  elliptic  operator.  A  more  involved  argument  shows  that 
the  fundamental  solution  satisfies  an  estimate  of  the  type  (2.8)  with  the  expression  x  +  at 
replaced  appropriately  by  solutions  of  ^  =  A,(x)  x(0)  =  x,  i  =  I,. . .  and  with  Cm,p 

possibly  growing  in  time  like  Cm.pfi*'*  for  k  fixed. 

Our  numerical  procedure  involves  the  compression  of  the  matrix  A",  which  for  the  purpose 
of  analysis  only,  we  shall  view  as  the  discretization  of  the  fundamental  solution  for  either 
(2.5)  or  (2.7), 

iA-)ik  =  G{xj,y,,n 
where  the  interval  [0, 1]  is  discretized  via 

xj  =  ^,  j  =  N  =  2\ 

[0, 1]  X  [0, 1]  is  discretized  via  {xj,yk),  and  t”  =  nAt  =  nAAx,  n  =  0, 1, . . .  . 

We  now  adapt  the  terminology,  notation,  and  results  of  [1]  to  this  unsteady  problem 

(1.1). 

Finite  difference  schemes  approximating  (1.1),  e.g.  (2.4)  are  regarded  as  acting  on  a 
vector  which  is  to  be  viewed  as  approximating  u(x,0)  on  the  finest  scale: 

s°  =  2^  J  ip(2''x  —  k  +  l)u(x,0)dx 

=  J  f{x)ip^k{x)dx. 

All  functions,  both  continuous  and  discrete,  are  extended  periodically: 

u{x,t)  =  u(x  +  l,t) 


ut  +  A{x)ux  =  C{x)u  +  (— l)^’^^(Ax)^'’  ^x 


,0  =  „o 

^k+N  — 


etc. 


The  function  <p  satisfies 

2m— 1 

~  P) 

p=0 

The  function  V’(a?)  which  will  generate  an  orthonormal  basis  is  obtained  via 

2m-l 

H  9p+M2x  ~  p) 

p=0 
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with  Qp  =  (—1)^  P  =  Ij  ■  •  •  j  2m  and  /  ifi{x)dx  =  1. 

The  coefficients  {/ip}p=i  are  generally  chosen  so  that 

'ipj,k{x)  =  2~^(2“-’x  -  k  +  l), 

for  j,  k  integers,  form  an  orthonormal  basis  and  in  addition,  the  function  ^(x)  has  m  van¬ 
ishing  moments 

j il>{x)x^dx  =  Q,  ^  =  0, 1, . . .  ,m  —  1. 

Also  we  define 

(fijk  =  2~3(p{2~^x  -  A:  -h  1). 

Finally,  we  assume  that  there  exists  a  real  constant  Tm{Ti  =  1)  following 

conditions  are  satisfied: 

J  (fi(x  +  Tm)x^dx  =  0  for  ^=l,...,m— 1, 

and  /  <p{x)dx  =  1. 

In  this  case  the  quadrature  formula  becomes: 

+  ocjv-'”*’))) 

and  the  initial  discretization  error  is  up  to  uniform  translation. 

The  decomposition  of  the  vector  {5°, . . .  ,S2i'}  into  the  basis  we  use  to  compute  with 
comes  via 

K}—  {4} -^{4}  - 

\  {4}\  •••\  {4}. 

This  is  implemented  in  0{N)  operations  using: 

p=2m 

4  =  K4+2k-l 

p=l 

p=2m 

4  =  E  gp^+lk-i 

p=i 

and  the  d\  are  viewed  as  periodic  sequences  with  period  2"“^. 

The  orthonormal  basis  consists  of 


[d\, . . .  ,d}iv,  d], . .  .d\, . . .  s”]. 


The  inverse  mapping  can  also  be  done  in  0{N)  operations. 
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Each  of  the  is  thought  of  as  approximating 

4  =  /  f{x)(pjkix)dx  = 

2-r-^)[f{2~^+^{k-l  +  Tm)) 

while  each  is  thought  of  as  approximating 

4  =  /  f{x)i>jk{x)dx. 

The  numerical  procedure  effectively  transforms  the  approximate  discretization  of  the 
matrix  G{xj,yk,t’')  which  is  {A^)jk.  Estimate  (2.8)  (corresponding  to  (4.5)  and  (4.6)  of  [1], 
uniform  in  all  parameters,  indicates  (via  an  argument  of  [1])  that  truncating  A"  by  removing 
elements  of  a  band  of  width  b  >  2m  around  a  shifted  diagonal  (and  its  periodic  extension) 
i.e.,  those  for  which 

\j  —  k  —  aAn|  >  b>  2m, 
which  replaces  A'^  by  A”'*,  leads  to  an  estimate 

for  C  depending  only  on  G. 

It  also  follows  easily  that  for  large  N  and  fixed  precision  e,  only  O(A^logA^)  elements 
will  be  greater  than  e.  Alternatively,  by  discarding  all  elements  that  are  smaller  than  a  fixed 
threshhold  we  compress  it  to  0{N log  N)  elements.  Again  following  the  discussion  in  [1],  we 
note  that  this  naive  approach  is  to  construct  the  full  matrix  in  the  wavelet  basis  and  then 
to  threshhold.  Clearly  this  is  an  0{N^)  operation. 

Since  we  have,  d  priori,  the  structure  of  the  singularities  of  the  matrix  A"  the  relevant 
coefficients  can  be  evaluated  by  using  the  quadrature  formulas.  Estimate  (2.8)  guarantees 
that  this  procedure  requires  0{N  log  N)  operations. 


Remark  R2.  It  is  interesting  to  note  that  so  called  unstable  difference  schemes  can  be  used 
without  any  drastic  loss  of  efficiency.  If  (2.1)  is  approximated  by, 

uf  ^  0/2, 


(2.9) 


vPj  =  uo{xj),  ;  =  l,2,...,iV 
the  algorithm  is  not  stable  for  any  fixed  A  >  0,  see  e.g.  [7]. 
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The  approximation  does  converge  if  At  <  CAx'^,  (A  <  CAx)  with  an  amplification 
factor  1  +  0(At).  The  number  of  timesteps  for  i  =  0(1)  calculation  will  be  large,  n  = 
0(Ai~^)  =  0{N^).  This  is  devastating  for  the  standard  explicit  algorithm  (1.2)  but  will 
only  aifect  the  complexity  of  (1.4)  by  a  constant  factor.  The  number  of  iterations  (m  in 
(1.4))  will  increase  from  log(iV)  to  log(iV^). 

Our  approach  is  in  general  not  as  favorable  for  multidimensional  hyperbolic  systems, 

a 

9tU  +  ^Ajix)dj:jU  =  f{x),  X  e  R‘^, 

(2.10) 

u(x,0)  =  uo(a^)- 

When  u  is  a  scalar  or  if  the  system  can  be  diagonalized  the  algorithm  (1.4)  works  well.  The 
solution  is  given  by  integration  along  characteristics  and  the  support  of  the  Green’s  function 
is  a  small  number  of  points  (see  Remark  (Rl)  above).  In  the  idealized  case  each  row  of 
A''  consists  of  a  fixed  number  of  6-functlons.  Its  wavelet  representation  will  have  log(A^‘^) 
nonzero  terms.  The  overall  complexity  for  (1.4)  is  then  0((log  A^)^A^‘^)  when  the  knowledge 
of  the  location  of  the  ^-functions  is  used.  This  is  better  than  the  standard  estimate. 

In  general,  however,  the  Green’s  function  for  (2.6)  has  a  support  with  positive  volume 
in  and  with  a  singular  support  of  positive  measure  in  Hausdorff  dimension  d  —  1.  The 
representation  of  the  singular  support  consists  of  C?(A/'‘^“^)d-functions  in  each  row  of  A*'. 
This  corresponds  to  0{\og{N)N'^~^)  wavelets  and  the  overall  algorithm  contains  at  least 
(C7(logiV)^iV^‘*"*)  wavelets. 

For  general  multidimensional  problems  the  new  algorithm  is  still  of  interest  in  special 
czises,  e.g.,  if  the  solution  is  needed  only  at  a  fixed  number  of  points  and  if  it  is  needed  for  a 
large  number  of  different  data  uq,  f. 


3.  Parabolic  Problems 


The  Green’s  function  for  parabolic  problems  is  smooth  in  contreist  to  the  hyperbolic  case. 
The  pure  initial  value  problem  for  the  heat  equation. 


(3.1) 


dtu  =  Au,  t  >  0,  X  €  R'^1 
u(x,0)  =  uo(x). 


hM  a  solution  of  the  form, 

(3.2)  u{x,t)  =  (47r<)"‘'/^  [  exp(-lx  -  y\^ /At)uo{y)dy. 

JR* 
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In  bounded  domains  the  kernel  has  to  be  changed  slightly  depending  on  the  boundary 
conditions.  For  positive  t(=  nAt)  each  row  in  is  always  an  approximation  of  segments  of 
regular  functions. 

Our  new  technique  is  in  general  more  favorable  for  parabolic  problems  than  hyperbolic 
ones.  The  structure  of  the  matrix  B  in  (1.4)  is  siiripler.  When  t  increases  the  kernel  becomes 
smoother  and  ajk  can  be  truncated  to  zero  for  all  k  when  j  is  lajge  enough. 

Explicit  methods  for  (3.1)  also  requires  more  operations  than  for  hyperbolic  problems 
when  the  standard  method  is  used.  This  follows  from  the  parabolic  stability  requirement, 

(3.3)  At  <  const.  lAarp. 

The  new  technique  is  only  marginally  affected  by  the  constraint  (3.3).  Compare  here  the 
discussion  above  for  unstable  hyperbolic  methods. 

In  more  general  higher  order  multidimensional  parabolic  cases  the  fundamental  solution 
of,  e.g., 

Ut  +  (— A)‘^u  =  0 
is 

I 

Gd{x,  0  =  ^  exp(i^  •  X  - 

This  is  merely  a  multidimensional  and  rescaled  version  of  the  fundamental  solution  used  in 
(2.8),  and  a  simpler,  but  multidimensional  version  of  (2.8)  is  just: 

\xr+'^D^Gd{x,t)\<Cn,d. 

Moreover  Cmd  is  arbitrarily  small  if  t  is  large  enough  (this  of  course  requires  the  nonexistence 
or  other  special  behavior  of  lower  order  terms). 

The  matrix  compression  technique  is  easy  here  (for  periodic  problems  without  boundary 
conditions)  because  the  significant  terms  of  [A*']  lie  near  the  main  diagonal  and  its  periodic 
extension  in  one  dimension.  In  two  space  dimensions  (as  is  usual  for  elliptic  operators),  we 
also  need  to  consider  diagonals  i  =  j  ±  kN  for  0  <  fc  <  d.  Recall  A  is  an  x  matrix  in 
2D. 

It  is  clear  that  a  priori  thresholding  (to  obtain  0(e)  precision)  near  the  image  of  these 
diagonals  will  give  us  an  0{N‘^{\ogN)^)  operation  for  each  evaluation  of  the  solution,  where 
d  is  the  number  of  space  dimensions  for  the  problem. 

4.  Numerical  Experiments 

The  algorithm  (1.4)  was  applied  to  hyperbolic  problems  in  one  space  dimensions  and  to 
one  and  two  dimensional  parabolic  problems.  Various  difference  approximations  and  wavelet 
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spaces  were  used.  We  present  results  concerning  the  accuracy  of  the  calculations  and  the 
sparsity  of 


4.1  Hyperbolic  problems.  Consider  the  following  scalar  hyperbolic  problem; 

dtu  +  a(x)dxU  =  f{x) 

(4.1a) 

u(x,0)  =  uo{x) 

with  periodic  boundary  conditions  (0  <  x  <  1).  We  made  the  following  choices; 
(4.1b)  a(x)  =  0.5  +  0.115  sin(47rx) 


(4,1c) 


/(x)  =  cos(47rx) 


(4. Id)  no(a:)  =  sin(4;rx). 

In  the  discretization,  Ax  =  1/1024  and  At /Ax  =  1.  The  wavelet  transform  operator 
S  uses  the  Daubechies-8  wavelets,  which  have  8  coefficients  and  have  4  vanishing  moments. 
Finite  difference  schemes  of  order  1, 2,3,4,  and  5  of  accuracy  are  tested. 

These  finite  diflference  schemes  are  obtained  as  follows.  In  each  interval 


(4.2)  /j,_i  =  {x/(i/  —  l)Ax  <  X  <  i/Ax} 

a  polynomial  of  degree  k  is  constructed.  This  polynomial  interpolates  the  two  points 
(x„_i,u"_j)  and  (x*,,it”)  and  A:  —  1  of  its  neighbors.  If  k  is  even  these  interpolation  points 
go  from  x^_*  to  If  k  is  odd  they  go  from  to  This  gives  us  a 

leconstruction  function  which  is  a  polynomial  of  degree  k  in  each  /^_i  and  is  continuous, 
but  generally  not  differentiable  at  the  boundary  points  x^-i  and  x^.  We  call  this  function 

ir*'*(x) 

To  approximate  (4,1)  at  the  grid  points  (x„,t”'''^)  we  solve  (4.1)  “exactly”  with  initial 
data 

(4.3)  tXA^(x,r)  =  i2”-*(x) 

for  t”  <  t  <  evaluate  the  solution  at  (x„,t”'''*),  and  set  =  nAi(x,/,t"‘''^).  We 
require  ME22^lEll<i^  so  the  solution  depends  only  on  data  in  I^_i  if  a(x)  >  0  and  if 
a(x)  <  0. 
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In  the  special  case  when  a{x)  =  a,  constant,  then 

yn+l  = 

-tn+l 

+  /  f{Xu  -  -  s 

Jt” 

In  the  case  when  /  =  0  we  get  some  familiar  schemes:  For  A:  =  1  this  is  just  the  first  order 
accurate  upwind  difference  scheme  (2.4).  For  A:  =  2  this  is  just  the  classical  Lax-Wendroff 
second  order  accurate  three  point  scheme,  see  e.g.  [7].  For  k  =  3,4,5  the  schemes  axe  less 
studied,  but  are  known  to  be  stable,  see  e.g.  [9]  and  the  references  therein. 

For  variable  coefficients  the  result  is 


(4.5a) 

where  x^{t)  solves 
(4.5b) 


(4.5c)  x,.(r+^)  =  X,,. 

A  fourth  order  Runge-Kutta  method  is  used  to  integrate  the  O.D.E.  (4i5b,c)  and  Simpson’s 
rule  is  used  to  evaluate  the  integral  in  (4.5a).  The  result  of  this  approximation  to  the  right 
side  of  (4.5a)  is  defined  to  be  u”’*"' . 

Returning  to  the  present  case  the  computations  ran  13  steps  until  t  =  4,  that  is, 
(5i45“^)^*’  was  computed. 

At  each  step  n  the  number  of  elements  of  A”  and  (jSAjS"^)’’  whose  absolute  values  are 
greater  than  10“'*  is  shown  in  table  1.  This  is  for  methods  whose  order  of  accuracies  go  from 
one  through  five.  The  results  are  also  plotted  on  Figure  1. 

These  significant  elements  are  located  near  the  sub-diagonal  corresponding  to  the  char¬ 
acteristic  curve  which  is  known  a  priori.  The  image  of  these  locations  in  shown 

on  figure  2,  has  total  length  of  0{N  log  N)  elements  where  N  =  1024. 

In  the  computation  of  (5A5“^)",  first,  from  the  knowledge  of  the  PDF,  we  figure  out  the 
structure  of  the  singularities  of  A  and  its  image  in  {SAS~^)".  Then  we  compute  (SAS"*)^”  = 
(5'A5'“*)"  ♦  (5'A5'"*)'*  considering  only  the  elements  in  a  neighborhood  of  the  singularities. 
In  particular,  we  define  the  neighborhood  of  a  singularity  to  be  locations  whose  distance 
from  the  singularity  are  less  than  or  equal  to  5.  If  the  singularities  lie  on  a  subdiagonal  and 
its  periodic  extension  its  neighborhood  form  a  subband  of  bandwidth  11  (the  wavelet  filters 
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have  8  elements).  This  bandwidth  is  independent  of  the  time  t  (the  step  n)  and  the  size  of 
the  problem.  The  errors  due  to  the  subband  truncation,  measured  by  ||u"  —  {t"ll/||u"||,  are 
shown  in  table  2b.  Table  2a  shows  the  relative  error  between  the  subband  truncation  and 
the  exact  solution.  Here  and  throughout,  “||  •  ||”  denotes  the  norm.  Table  2c  shows  the 
relative  error  between  the  subband  truncation  and  untruncated  under  grid  refinement  for  the 
various  orders.  Unsurprisingly,  since  the  relative  length  of  the  subband  which  is  preserved 
decreases  linearly  with  grid  size,  the  error  increases,  but  only  slightly  under  this  process. 

We  note  that  the  compression  (as  seen  in  Figure  1  and  Table  1)  is  better  for  odd  order 
than  for  even  order  schemes.  This  is  perhaps  not  surprising  since  (2.7)  models  schemes  of 
odd  order  accuracy.  Singularities  behave  a  bit  differently  for  even  order  (say  order  =  2p) 
schemes.  These  are  modeled  by 


(4.6) 


ut  +  aux 


=  4(Ax)2p 


u 


+i-lYkpiAxf^+^ 


u 


where  kp  >  0  and  £p  are  nonzero  constants.  The  odd  order  dispersive  term  above  may  tend 
to  spread  singularities  of  the  fundamental  solution  spuriously. 

Finally  table  3  shows  the  relative  error  due  to  truncation  when  the  band  width  of  the 
subband  is  9,  11,  and  13  for  the  methods  of  first  and  second  order.  Figures  3a  and  3b 
compare  the  truncated  versus  the  approximate  solutions  due  to  truncation  of  bandwidth  9 
for  the  first  and  second  order  methods  (the  truncated  graphs  are  dotted). 


4.2  Unstable  Schemes.  For  theoretical  interest,  we  apply  the  method  to  a  finite  difference 
scheme  which  is  unstable  for  ^  =  A  >  0 

(4.7a)  uf'  =  <  -  -  uU)/2, 


(4.7b)  uj  =  u<i(ai). 

The  amplification  factor  of  this  scheme  is 
(4.8)  1  —  Ai  sin0  =  r(e’^),  — tt  <  0  <  1 

so 

lr(e’®)|  =  (l  +  A^sin^e)*. 


At  <  2c{AxY 


This  means  that  if 
(4.9) 
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for  some  c  >  0,  then 

(4.10)  |1A"11p  < 

The  restriction  (4.9)  means  that  the  operation  count  for  this  explicit  method  would  be 
O(N^)  if  we  were  silly  enough  to  use  it.  However  our  compression  method  allows  for  an 
operation  count  of  0(iV(log  N)^)  for  the  reasons  described  above. 

Table  4  shows  the  number  of  elements  in  A”  and  (SAS~^)'*  whose  absolute  values  are 
greater  than  10“^.  We  choose  a  bigger  threshold  here  since  we  took  =  1  and  nAt  =  2, 
so  ||A”||,  as  estimated  in  (4.10)  grows  to  be  roughly  10  when  we  are  finished  computing. 

The  error  as  measured  by  (subband  truncation  using  bandwidth  11)  was  0.0136. 

We  also  performed  convergence  studies  as  we  refined  the  grid  for  this  method.  Figures 
(4a,b,c)  compare  the  numerical  (untruncated)  using  dots  versus  exact  solution  for  m  = 
128, 256, 512  grid  points.  The  result  indicates  a  second  order  method,  as  it  should,  since  At  = 
(Ax)*.  Figures  (5a,b,c)  compare  the  truncated  bandwidth  (using  dots)  vs  the  untruncated 
for  this  method  for  m  ^  128,256,  and  512  grid  points. 

The  relative  error  decreases  with  mesh  refinement.  The  truncation  error  equation  associ¬ 
ated  with  this  scheme  involves  limited  antidiffusion.  Perhaps  this  accounts  for  this  behavior. 


4.3  System  of  Hyperbolic  Equations.  We  apply  the  method  to  solving  the  system  of 
hyperbolic  equations: 


(4.11o) 


u  a  0  « 

+  n 

w  0  —a 


on0<x<l,t>0  with  the  boundary  conditions  and  initial  conditions: 

v(0yt)  =  w{0,  t) 

w{l,t)  =  v{l,t) 

(4.11b) 

u(x,0)  =  uo(x) 

tu(x,0)  =  iwo(a^) 

the  coefficient  a  is  chosen  to  be  constant: 


a  =  0.115. 


The  numerical  method  used  is  the  first  order  accurate  upwind  method  described  above. 
The  results  are  similar  to  the  scalar  case,  except  the  structure  of  the  singularities  in  the 
matrices  is  more  complicated.  We  have  to  keep  track  of  reflections  of  singularities  at  the 
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boundaries  which  is  quite  simple  in  this  case.  The  number  of  elements  in  A"'  and 
whose  absolute  values  are  greater  than  10“'*  is  shown  on  table  5,  and  is  plotted  on  figure  6. 
The  relative  error  due  to  the  subband  of  width  11  truncation,  measured  by  Hu"  —  {t"|j/||u"||, 
is  0.0149. 

The  structure  of  the  elements  whose  absolute  values  are  greater  than  10“*  of  and 
is  shown  in  figures  (7a, c),  while  Figure  {7b)  shows  the  image  of  a  subband  of 
bandwidth  11  in 


4.4  Parabolic  Problems.  We  do  experiments  on  the  following  parabolic  problem: 

dtu  =  dx{a(x)dxu)  +  f{x) 

(4.12) 


u(x,0)  =Uq{x) 

with  periodic  boundary  conditions  (0  <  x  <  1).  We  made  the  following  choices: 
a{x)  =  0.5  +  0.25sin(2s'x) 

/  (x)  =  —  TT^  cos(27rx)^  +  7r^(0.5  +  0.25  sin(27rx))  sin(27rx) 


uo((e)  =  sin(47rx). 


The  discrete  setting  and  the  wavelets  are  the  same  as  in  the  hyperbolic  problem.  We  use 
the  simple  explicit  central  difference  scheme  (4.13) 

(4.13) 

+ 

where 

with  Al/(Ax)^  =  0.25.  The  number  of  significant  elements  in  and  {SAS~^)'^  is  shown  on 
table  6,  and  is  plotted  on  figure  8. 

For  the  parabolic  problem,  the  large  elements  of  A  are  in  the  neighborhood  of  the  main 
diagonal.  Their  wavelet  transform  image  is  shown  in  figure  9.  The  relative  error  due  to 
subband  truncation  weis  0.0025. 


4.5  Two-dimensional  Parabolic  Problems.  We  consider  the  following  problem: 

dtU  =  O.iidxx'^  +  *io,\2dxyU  +  Cl229yyU 
u{x,y,Q)  =  uo{x,y) 
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with  periodic  boundary  conditions  (0  <  i  <  1,  0  <  y  <  1).  We  choose 

^nix,y)  =  0.5  +  0.25  8in(2irx) 

<^\2{x,y)  =  0.115  sin(27rx)cos(2fl-y) 

<^22(x,  y)  =  0.5  +  0.25  co3(2iry) 

^o(a!,y)  =  sin(47rx)  +  cos(87rx). 

We  use  a  standard  two-dimensioncil  explicit  central  dilFerence  scheme.  The  two-dimensional 
data  j  =  k  =  I...N2  forms  a  one-dimensional  vector  in  the  following  way 

{ui.i  . .  .  Ui.iVj,  U2,l  . . .  «2,Na,  •  •  • ,  . . .  USi.Nt}- 

To  reduce  the  size  of  the  problem,  N2  is  much  less  than  Wi.  In  particular  we  took  Ni  = 
128,  7^2  =  8  that  is.  Ax  =  Ay  =  |. 

The  compression  worked  quite  well.  Table  7  shows  the  number  of  elements  in  A”  on 
(5A5“^)"  whose  absolute  values  are  greater  than  10~^.  The  relative  error  due  to  subband 
truncation  was  0.0066. 
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Figur<  3*  Thincation  vmus  oontruncatad  apprudmate  solution,  first  order  method  trun¬ 
cated  at  bandwidth  9  (Truncated  is  dotted^ 
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Figure  3b;  Truncated  versus  nootruncated  ^preadmale  aolution,  accond  order  method, 
truncated  at  bandwidth  9  (Truncated  is  dotted). 
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Figont;  SyiUm  cf  hypcbolit  «qu«tioB»'.  Tht  aaaibtr  rf  ■ImiwiU  in  «Bd  {SA5~*  )* 
wfagit  abMiuit  laiiMi  aic  pottr  thia  t0~^. 


Figun  7«:  Sjmcai  of  kypcrbolic  •quitiooi  pattans  of  ■gniBniit  fa 

A',  It  ■  2046 


rtfortTb:  >yiwrflorp»fcalic«gMti— p»tf»rfiignttc»ntii«D»<far(545~*)*t "  ■ 

2041.  iMgt  af  bMdvidtb  >1  HtMDd  riagaitf  ngp«t 
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